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Abstract 

Quantum magnetism in low dimensions has been one of the central areas 
of theoretical research for many decades now. One of the key reasons for the 
long standing interest in this field has been the existence of simplified models, 
which serve as paradigms for understanding the role of strong interactions 
in many-electron systems. Although simple, these models quite often can 
not be solved exactly. In this review, we discuss a variety of analytical and 
numerical methods, which treat the system in a systematic and controlled 
manner. The central method employed in all the studies is the density matrix 
renormalization group (DMRG) method. This is supported by small scale 
numerical exact calculations, and by analytical methods using field theoretic 
techniques. We have considered a number of magnetic systems, from magnetic 
clusters to extended lattices, and have found some novel quantum ground 
states and low-energy elementary excitations. In some cases, we have also 
employed the finite-temperature DMRG method to accurately compute the 
low-temperature thermodynamic properties such as specific heat and magnetic 
susceptibility. 
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1 Introduction 



The nonrelativistic Schrodinger equation of a system of electrons is spin independent. 
It therefore appears at first glance that the solutions of the Schrodinger equation 
should also be spin independent. However, the indistinguishability of the electrons 
forces the total wave function which is a product of the spin wave function and 
the spatial wave function to be antisymmetric. This in turn implies that for two 
electrons, a spatially symmetric wave function should be associated with an anti- 
symmetric spin wave function and vice versa. The different charge distribution in 
the spatially symmetric and antisymmetric wave functions leads to different coulomb 
repulsions by virtue of which the spin states which are symmetric and antisymmetric 
have different energies fl]]. Dirac represented the splitting between the energies of 
the two spin states by the spin operator — 2 JSj ■ Sj, where J is the exchange integral 
involving the two spatial orbitals in which the two electrons are singly occupied. In 
most open shell atomic systems the exchange integral J is large enough to force the 
total spin of the ground state configuration to be the largest permissible value. This 
in essence is the Hund's rule of maximum multiplicity and is also the reason why 
we find transition and rare earth metal ions in the high spin states in nature. 

In solids containing transition metal or rare earth ions surrounded by ligands, 
the relative alignment of the unpaired spins at the metal site is not at all obvious. 
In order to understand this, we should examine the possible pathways for the de- 
localization of the valence electrons in the system. If the favorable derealization 
pathways involve antiparallel alignment of the metal ion spins then the nature of 
the exchange interaction between the metal ion spins is antiferromagnetic and fer- 
romagnetic otherwise. This follows from the fact that derealization of electrons 
lowers kinetic energy of the electrons and therefore the ground state corresponds to 
an alignment of spins that allows for maximum derealization. This is indeed the 
reason why the ground state of a hydrogen molecule is a spin singlet. In a system 
with degenerate partially occupied orbitals, the Hund's coupling favors high-spin 
alignment of the electrons on an ion. If derealization pathways exist that allow 
for these high-spin states in the process of derealization, then the alignment of 
the spins on two neighboring centers will be ferromagnetic. If derealization path- 
ways exist only when these high-spin states are aligned antiparallel, one would have 
an antiferromagnetic alignment of the spin. Thus, the overall nature of the spin 
alignment is governed by a competition between the Hund's coupling and electron 
derealization |fj, [|. 

Given a collection of spins, the exchange Hamiltonian for the system is written 
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as 

H = J ijSi ' S i > (!) 

ij 

where Jjj is the effective exchange integral for the interaction between the spins at 
sites 'i' and 'j'. Since the spins in the cluster arise from unpaired electrons of a 
transition metal ion in a crystal field, it is natural to expect that the spin-orbit as 
well as spin-spin interactions of the electrons in the ion could alter the nature of the 
total spin on the ion by giving the net spin a preferred direction of orientation. Such 
a situation can be easily handled by treating each J^- as a vector and generalizing 
the exchange Hamiltonian as 

TT _ / JX qX qX i TV OP OV I JZ qz qz \ /r>\ 

n — V J ij D i D j i J ij>~>iOj T J i jO i Dj ) . [Z) 

Such a model is often referred to as the XYZ spin model . Two extreme cases are 
often studied, (i) the spin is assumed to have no projection on the X-Y plane, in 
which case the resulting model is the Ising model and corresponds to scalar spins, 
and (ii) the spin is assumed to have no projection on the z-axis, in which case we 
have an XY model or a planar spin model. The Ising model is a discrete classical 
model as it consists of no noncommuting operators in its Hamiltonian while the XY 
model could be classical or quantum mechanical. Usually, while dealing with large 
site spin systems, it is not uncommon to assume that the spins are classical, in the 
spirit of Bohr's correspondence principle. 

In the crystalline state, the spins in the solid would be arranged on a lattice. 
If the exchange interaction is predominant between spins along a single crystalline 
direction, the model could be treated as a one dimensional array of spins. There are 
many examples of solids in which this holds [||. Likewise, it is also possible that the 
interactions amongst spins is large along two crystallographic directions and weak 
along a third direction and this would result in two-dimensional spin system |J . 

In this review article, we will mainly concern ourselves with the study of isotropic 
spin clusters and one-dimensional spin systems, sometimes in the presence of an ex- 
ternal magnetic field. We will be mostly interested in properties of the ground 
state and the low-lying excitations, since these are the states which govern the 
low-temperature properties of systems such as the specific heat and magnetic sus- 
ceptibility. The symmetries of a system often enable us to characterize the energy 
eigenstates in terms of quantum numbers such as the total spin S to t, the component 
of the total spin along some particular direction, say, S to t,zi s P m parity (which is a 
symmetry for states with S tot ^ z = 0), the wave number k for a translation invariant 
system, and possibly other spatial symmetries depending on the structure. We will 
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discuss below how the use of symmetries can help to lessen the numerical effort 
required to study the low-energy states. 

In the next two sections we introduce some exact numerical methods, and de- 
scribe the application of these methods to magnetic clusters. In section 4, we discuss 
two analytical methods which use field theoretic approximations. In section 5, we 
describe an innovative way of solving the spin Hamiltonians, which goes beyond the 
conventional techniques and is based on the density matrix renormalization group 
(DMRG) theory. Various applications of DMRG to the properties of low-dimensional 
extended chains are described in sections 6 to 8. 

2 Exact Calculations 

The properties of a spin Hamiltonian can be computed from the eigenstates of the 
Hamiltonian which are in turn obtained by setting up the Hamiltonian matrix in 
a suitable basis and diagonalizing it thereafter. While the procedure itself is quite 
straightforward, the space spanned by the Hamiltonian rapidly increases with the 
number of the spins in the system. The Fock space dimensionality of a system of n 
spins with spin Sj is given by 

D F = f[ (2s t + I) . (3) 

The Hamiltonian matrix is block-diagonal in structure with each block corresponding 
to specified values of the quantities conserved by the Hamiltonian. Thus, for an 
isotropic spin system, the z-component of the total spin, Ms, and the total spin S are 
conserved. Restricting the Fock space to specified values of M s and S gives Hilbert 
spaces whose dimensionalities are smaller than the Fock space dimensionality. 

While constructing spin basis functions which are eigenstates of the total S z 
operator is quite simple, construction of spin adapted functions (SAF, eigenstates of 
the total S operator) is not direct. Perhaps the simplest and chemically appealing 
way of constructing SAF is by the valence bond (VB) method which uses the Rumer- 
Pauling rules. This method is best illustrated by applying it to a system of 2n spins, 
each possessing a spin of half, in the total spin S = sector. A total spin singlet 
can be formed by choosing pairs of sites and spin coupling each of them to obtain a 
singlet. The product of these singlet pairs will be a spin eigenstate of the operator 
^totai- This is illustrated in Fig. 1. However, there are more ways to spin couple 
in pairs than the number of linearly independent singlet states, e.g., the state |3 > 
in Fig. 1 can be expressed as a linear combination of the states |1 > and \2 >. 
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The overcompleteness can be avoided by resorting to the Rumer-Pauling rules. To 
implement this rule, we arrange the 2n spins at the vertices of a regular 2n-gon. We 
draw lines between pairs of sites that are singlet paired. According to the Rumer- 
Pauling rule, the subset of these, encompassing all diagrams (to be called 'legal' 
diagrams) with no crossing lines, forms a complete and linearly independent set of 
states |J. 

The Rumer-Pauling rules can be easily extended to construct complete and lin- 
early independent basis sets in higher spin Hilbert spaces involving spin- 1/2 objects. 
This is done most easily with the help of phantom sites. If we wish to construct VB 
diagrams for total spin S subspace involving n spin- 1/2 objects, then we introduce 
2S additional sites to be called phantom sites. Besides imposing the Rumer-Pauling 
rules on the diagrams with n + 2S sites, we impose the additional constraint that 
there should be no singlet lines amongst the 2S phantom sites. In Fig. 2, we show 
a few examples of VB diagrams with higher total spin. 

It is also quite simple to extend the VB rules to spin clusters made up of different 
site spins M. If the spin at a site is Sj, then we replace this site by a set of 2si sites, 
each with spin-1/2. We then proceed with constructing the VB basis, as though 
the system is made up entirely of spin-1/2 objects, with one difference, namely, we 
impose the additional constraint that there should be no singlet lines within the 
subset of 2si sites which replace the spin Sj at site %. The VB diagrams with total 
spin S > | are constructed as before with the help of phantom sites. An example 
of a legal VB diagram involving higher site spins is shown in Fig. 2. 

Generating and storing the VB diagrams on a computer is also quite simple. 
We associate one bit with every site. The state of the bit is '1', if in the VB 
diagram a line begins at the site and the state is '0' if a line ends at the site. Thus, 
we can associate an integer of n bits with every VB diagram involving n spin-1/2 
objects. This association is unique if we decipher the bit pattern of the integer 
corresponding to the diagram from inside-out, much like expanding an algebraic 
expression with multiple parentheses. In Fig. 2, we have also shown the bit pattern 
and the associated integer for each VB diagram. The VB diagrams are generated on 
a computer by checking the bit pattern in all n-bit integers to see if they satisfy the 
criterion for representing the desired VB diagram. This also allows us to generate 
the VB diagrams as an ordered sequence of the integers that represent them, a fact 
that helps in rapid generation of the Hamiltonian matrix. 

The Hamiltonian matrix in the VB basis can be easily constructed by knowing 
the action of the operator Sj-Sj for spin-1/2 particles (i) on a singlet line joining sites 
% and j and (ii) on the pair of sites % and j singlet paired to two different sites %' and f 
(Fig. 3). In the case of Hilbert spaces with nonzero total S, the Hamiltonian involves 
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spin exchange between the real sites only. However, these exchange operators could 
lead to VB diagrams in which the phantom sites are interconnected. In this event, 
simply neglecting these resultant states is sufficient to ensure that we are dealing 
exactly with the spin S Hilbert space. The Hamiltonian for spin clusters with 
arbitrary spins can be treated as consisting of operators with only spin-1/2 objects. 
This is done by replacing the spin exchange operator between sites i and j, Sj • Sj by 
the operator (E^iTfc) ■ (S^Tj), where the operators f k and % are the usual spin-1/2 
operators. 

The matrix representing the Hamiltonian in the VB basis is in general nonsym- 
metric since the VB basis is nonorthogonal. However, the matrix itself is sparse. 
There exist efficient numerical algorithms || for obtaining the low-lying eigenstates 
of sparse nonsymmetric matrices, and it is possible to solve a nonsymmetric matrix 
eigenvalue problem for a million by million matrix with about a 100 million nonzero 
matrix elements on a powerful PC based workstation. 

While VB theory guarantees spin purity in the computed eigenstates of the spin 
conserving Hamiltonians, it has several drawbacks. Computation of quantities such 
as spin-spin correlation functions and spin densities are not easy because a site-spin 
operator operating on a VB diagram spoils the spin purity of the diagram. This 
could be overcome by converting the eigenstate of the Hamiltonian in the VB basis 
to the constant Ms basis. Another difficulty with the VB procedure is exploiting 
the spatial symmetries of the problem. Operation by a spatial symmetry operator 
on a legal VB diagram could lead to illegal VB diagrams. Disentangling these illegal 
diagrams into legal diagrams can be computationally prohibitive ||. 

In more general spin problems, it is often advantageous to use the constant Ms 
basis and exploit all the spatial symmetries. Partial spin symmetry adaptation in 
these cases is also possible by using the spin parity operator. The effect of the spin 
parity operator on a basis state is to flip all the spins in the state. In the Ms = sec- 
tor, it is possible to factor the Hilbert space into odd and even parity Hilbert spaces. 
The odd (even) parity Hilbert space is spanned by basis vectors with odd (even) 
total spin. This also has the effect of reducing the dimensionality of the Hilbert 
space besides providing partial spin symmetry adaptation. It is rather simple to 
set up the Hamiltonian matrix in the symmetry adapted basis. The Hamiltonian 
matrix is symmetric and usually very sparse. The lowest few eigenstates can be 
easily computed by employing the Davidson algorithm. Given these eigenstates, 
the computation of properties can proceed by converting an eigenstate in the sym- 
metrized basis into that in the unsymmetrized basis. The orthogonality of the basis 
states as well as the simple rules involved in obtaining the resultant when a basis 
state is operated upon by any type of spin operators in any combination affords easy 
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computation of a variety of properties of a magnetic system. 

The exact diagonalization techniques discussed above are in general applicable to 
systems whose Hilbert space dimensionality is about 10 million. The major problem 
with exact diagonalization methods is the exponential increase in dimensionality of 
the Hilbert space with the increase in the system size. Thus, the study of larger 
systems becomes not only CPU intensive but also memory intensive as the number of 
nonzero elements of the matrix increases rapidly with system size. With increasing 
power of the computers, slightly larger problems have been solved every few years. 
To illustrate this trend, we consider the case of the spin-1 Heisenberg chain. In 1973, 
ten years before the Haldane conjecture, De Neef 0] used the exact diagonalization 
procedure to solve a 8-site spin- 1/2 chain. In 1977, Blote 
Hamiltonian of a chain of 10 sites. In 1982, Botet and Jullien 
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diagonalized the 
increased this to 

12 sites. In 1984, Parkinson and Bonner jl3| solved the 14-site spin-1 problem. In 
the same year, Moreo [H] solved the 16-site spin-1 chain. In 1990, Takahashi [15 
pushed this up to 18 sites. And in 1994, Golinelli et al. JT^j have solved for the 
low- lying states of a 22-site spin-1 chain. The growth in chain length of the longest 
spin-1 chain solved is almost linear with time, roughly increasing by 2 sites in every 
three years. Just to remind ourselves, the Fock space dimensionality in this case 
increases as 3^ with chain length N. The size of the matrix also increases similarly 
and the CPU and storage scales quadratically with the size of the matrix, if we 
are targeting only a few eigenstates. For this reason, for systems which span much 
larger spaces, the focus has shifted to approximate techniques. 



3 Applications to Spin Clusters 

In recent years, some of the magnetic clusters that have been studied extensively are 
the Mn 12 []T7|, Fe 8 JH| and V15 |B| clusters. These clusters show many interesting 



phenomena such as quantum resonant tunneling and quantum interference [p0 
Basic to a proper understanding of these phenomena is a knowledge of the low- 
energy excitation spectrum in these systems. The methods discussed under exact 
diagonalization schemes allow us to calculate the low-energy excitation spectrum, 
given a set of exchange constants. However, the exchange constants themselves are 
not known with any certainty. Therefore, it is all the more important to be able to 
carry out exact diagonalization studies of low-lying states to infer the possible sign 
and magnitude of the exchange constants 0. 

For the Mni2 cluster, we show the geometry and the exchange parameters in 
Fig. 4. The crystal structure suggests that the exchange constant J\ is largest 



7 



and antiferromagnetic in nature ||21|| . Based on magnetic measurements, it has 
been suggested that J\ has a magnitude of 215K. The magnitude and sign of the 
other exchange constants are based on comparisons with manganese systems in 
smaller clusters. It has been suggested that the exchange constant J 2 and J3 are 
antiferromagnetic and have a magnitude of about 85K. However, for the exchange 
constant J4, there is no concrete estimate, either of the sign or of the magnitude. 
In an earlier study, the Mn m - Mn IV pair with the strongest antiferromagnetic 
exchange constant was replaced by a composite spin-1/2 object P3 and the exchange 
Hamiltonian of the cluster solved for three different sets of parameters. It was found 
that the ordering of the energy levels were very sensitive to the relative strengths 
of the exchange constants. In these studies, J4 was set to zero and the low-lying 
excited states were computed. Besides, only states with spin 5 up to ten could be 
obtained because of the replacement of the higher spin ion pairs by the composite 
spin-1/2 objects. 

The technique described earlier, however, allows an exact computation of the low- 
lying states of Mn 12 . The results of the exact calculations are presented in Table 
1. We note that none of the three sets of parameters studied using an effective 
Hamiltonian, gives the correct ground and excited states, when an exact calculation 
is performed. It appears that setting the exchange constant J 4 to zero, cannot yield 
an 5 = 10 ground state (Table 1, cases A, B and C). When J 3 is equal to or slightly 
larger than J 2 (cases A and B, Table 1), we find a singlet ground state, unlike the 
result of the effective Hamiltonian in which the ground state has S — 8 and 5 = 
respectively. The ground state has spin 5 = 6, when J3 is slightly smaller than J 2 
(case C, Table 1). In all these case, the first few low-lying states are found to lie 
within 20K of the ground state. 

When we use the parameters suggested by Chudnovsky |17| (case D, Table 1), we 
obtain an 5 = 10 ground state separated from an 5 = 9 first excited state by 223K. 
This is followed by another 5 = 9 excited state at 421K. Only when the exchange 
constant J4 is sufficiently strongly ferromagnetic (case E, Table 1), do we find an 
5 = 10 ground state with an 5 = 9 excited state separated from it by a gap of 35K, 



which is close to the experimental value |f23" l- The second higher excited state has 



5 = 8, and is separated from the ground state by 62K. 

In Fig. 5, we show the spin density [^| for the Mn 12 cluster in the ground state 



for the 5 = 10, Ms = 10 state. While the manganese ions connected by the strong 
antiferromagnetic exchange show opposite spin densities, it is worth noting that the 
total spin density on these two ions is 0.69, well away from a value of 0.5 expected 
if these ions were indeed to form a spin-1/2 object. 

The Fes cluster is shown in Fig. 6. Each of the Fe ions has a spin of 2 and 



S 



the ground state of the system has a total spin S = 10, with S — 9 excited state 
separated from it by about 20K. All the exchange interactions in this system are 
expected to be antiferromagnetic. While the structure of the complex dictates that 
the exchange interaction J2 along the back of the butterfly should be small in com- 
parison with the interaction J\ across the wing , in earlier studies it was reported 



that such a choice of interaction parameters would not provide a S = 10 ground 



state 26 



Results from the exact calculations of the eigenstates of the Fes cluster using 
three sets of parameters is shown in Table 2. In two of these cases, J2 is very much 
smaller than J\. We find that in all these cases, the ground state has a spin S = 10 
and the lowest excited state has spin, S — 9. One of the main difference we find 
amongst the three sets of parameters is in the energy gap to the lowest excited state 
(Table 2). For the set of parameters used in the earlier study, this gap is the lowest 
at 3.4K. For the parameter sets 1 and 3 f27|, this gap is respectively 13. IK and 



39. 6K. While in cases 1 and 2, the second excited state is an S = 8 state, in case 3, 
this state also has spin 9. 

The spin densities in all the three cases for the ground state are shown in Fig. 7. 
The spin densities in all cases are positive at the corners. In cases 1 and 2, the spin 
density on the Fe ion on the backbone is positive and negative on the remaining 
two Fe sites [pq| . However, in case 3, the negative and positive spin density sites 
for Fe ions in the middle of the edges is interchanged. This is perhaps due to the 
fact that in cases 1 and 2, the exchange constant J3 is less than J4, while in case 3, 
this is reversed. Thus, a spin density measurement can provide relative strengths of 
these two exchange constants. In all the three case, the difference between the spin 
densities in the ground and excited states is that the decrease in the spin density in 
the excited state is mainly confined to the corner Fe sites. 



4 Field Theoretic Studies of Spin Chains 

One-dimensional and quasi-one-dimensional quantum spin systems have been stud- 
ied extensively in recent years for several reasons. Many such systems have been 
realized experimentally, and a variety of theoretical techniques, both analytical and 
numerical, are available to study the relevant models. Due to large quantum fluc- 
tuations in low dimensions, such systems often have unusual properties such as a 
gap between the ground state and the excited states. The most famous example 
of this is the Haldane gap which was predicted theoretically in integer spin Heisen- 



berg antiferromagnetic chains 2G , and then observed experimentally in a spin-1 
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system Ni(C 2 H 8 N 2 ) 2 N0 2 (C104) ^0|. Other examples include the spin ladder sys- 
tems in which a small number of one-dimensional spin- 1/2 chains interact amongst 
each other It has been observed that if the number of chains is even, i.e., 



if each rung of the ladder (which is the unit cell for the system) contains an even 
number of spin-1/2 sites, then the system effectively behaves like an integer spin 
chain with a gap in the low-energy spectrum. Some two-chain ladders which show 
a gap are (VO) 2 P 2 7 |2|, SrCu 2 3 @ and Cu 2 (C 5 H 12 N 2 )2Cl4 0. Conversely, a 



three-chain ladder which effectively behaves like a half-odd-integer spin chain and 
does not exhibit a gap is Sr 2 Cu30s p3fl . A related observation is that the quasi-one- 
dimensional system CuGeOs spontaneously dimerizes below a spin-Peierls transition 



temperature then the unit cell contains two spin-1/2 sites and the system is 
gapped. 

The results for gaps quoted above are all in the absence of an external magnetic 



field. The situation becomes more interesting in the presence of a magnetic field |36 
Then it is possible for an integer spin chain to be gapless and a half-odd-integer spin 
chain to show a gap above the ground state for appropriate values of the field [37- 
45]. This has been demonstrated in several models using a variety of methods such 
as exact diagonalization of small systems and bosonization [|6 , f!7|| . In particular, 



it has been shown that the magnetization of the system can exhibit plateaus at 
certain nonzero values for some finite ranges of the magnetic field. Further, for a 
Hamiltonian which is invariant under translation by one unit cell, the value of the 
magnetization per unit cell is quantized to be a rational number at each plateau 



37J. In section 8, we will study the magnetization plateau which can occur in a 
three-chain ladder. 

In the next two subsections, we will discuss some field theoretic methods which 
can be used for studying spin chains and ladders. These methods rely on the idea 
that the low-energy and long- wavelength modes of a system (i.e., wavelengths much 
longer than the lattice spacing a if the system is defined on a lattice at the micro- 
scopic level) can often be described by a continuum field theory. 

4.1 Nonlinear a- model 

The nonlinear cx-model (NLSM) analysis of antiferromagnetic spin chains with the 
inclusion of J 2 (next-nearest neighbor coupling) and 5 (dimerization) proceeds as 



written as 



follows ||48|| . The Hamiltonian for the frustrated and dimerized spin chain can be 
H = E [ 1 ~ 5 ] ■ S m + J 2 £ S 4 ■ S l+2 . (4) 
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The interactions are schematically shown in Fig. 8. The region of interest is defined 
by J% > and < 5 < 1. We first do a classical analysis in the S — * oo limit 
to find the ground state configuration of the spins. Let us make the general ansatz 
that the ground state is a coplanar configuration of the spins with the energy per 
spin being equal to 



e 



S 2 



y (l + 5)cos0i + — (l-5)cosfl 2 + J 2 COS(0! 



6o 



(5) 



where 6\ is the angle between the spins S21 and S21+1 and 62 is the angle between 
the spins S2« and S2i_i. Minimization of the classical energy with respect to the 9i 
yields the following three phases. 

(i) Neel phase: This phase has Q\ = 6 2 = it; hence all the spins point along the 
same line and they go as 
1 -5 2 > U2/J1. 

(ii) Spiral phase: Here, the angles 6\ and 9 2 are given by 



TITI ' ' ' along the chain. This phase is stable for 



cos 81 
and cos 9 2 



1 



1 + 5 
1 

1-5 



1-5 2 



4J 2 



4/2M 1 + 5 2 J x 
1 - 5 2 5 4 J 2 



4J 2 /Ji 



1-5 2 Ji 



(6) 



where tc/2 < Q\ < 71 and < 8 2 < 6\. Thus the spins lie on a plane. This phase is 
stable for 1 - 5 2 < 4J 2 M < (1 - 5 2 )/5. 

(iii) Colinear phase: This phase (which needs both dimerization and frustration) is 
defined to have Q\ = tt and 62 = 0; hence all the spins again point along the same 
line and they go as • ■ • till • • • along the chain. It is stable for (1 — 5 2 )/5 < 4J 2 / J\- 

These phases along with their boundaries are depicted in Fig. 9. Thus even in the 
classical limit S — ► 00, the system has a rich ground state 'phase diagram' |49[| . 

We can now go to the next order in 1/S, and study the spin wave spectrum 
about the ground state in each of the phases. The main results are as follows. In 
the Neel phase, we find two zero modes, i.e., modes for which the energy Uk vanishes 
linearly at certain values of the momentum k, with the slope dujk/dk at those points 
(the velocity) being the same for the two modes. In the spiral phase, we have 
three zero modes, two with the same velocity describing out-of-plane fluctuations, 
and one with a higher velocity describing in-plane fluctuations. In the colinear 
phase, we get two zero modes with equal velocities just as in the Neel phase. The 
three phases also differ in the behavior of the spin-spin correlation function S(q) = 
X)„(S • S n ) exp(— iqn) in the classical limit. S(q) is peaked at q = {61 + 9 2 )/2, i.e., 
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at q = 7r in the Neel phase, at7r/2<g<7rin the spiral phase and at q = tc/2 in 
the colinear phase. 

To study the interactions between the spin waves, it is convenient to derive 
a semiclassical NLSM field theory which can describe the low-energy and long- 
wavelength excitations. The field theory in the Neel phase is given by a 0(3) NLSM 
with a topological term [^, |T7[ . The field variable is a unit vector <fi with the 
Lagrangian density 

where c = 2S(l-4J 2 -5 2 f' 2 is the spin wave velocity, g 2 = 2/[S(l-AJ 2 -5 2 ) 1/2 } is 
the coupling constant (which describes the strength of the interactions between the 
spin waves), and 9 = 2ttS(1 — 5) is the coefficient of the topological term (the integral 
of this term is an integer which defines the winding number of a field configuration 
4>(x, t)). Note that 9 is independent of J 2 in the NLSM. (Time and space derivatives 
are denoted by a dot and a prime respectively). For (#,mod27r) = 7r and g 2 less 
than a critical value, it is known that the system is gapless j|7|, 5T)f . For any other 



value of 9, the system is gapped. For J 2 = 5 = 0, one therefore expects that integer 
spin chains should have a gap while half-odd-integer spin chains should be gapless. 
This is known to be true even for small values of S like 1/2 (analytically) and 1 
(numerically) although the field theory is only derived for large S. In the presence 
of dimerization, one expects a gapless system at certain special values of S. For 
S = 1, the special value is predicted to be 5 C = 0.5. We see that the existence 
of a gapless point is correctly predicted by the NLSM. However, as we will see 
later, according to reliable numerical results from DMRG, 5 C is 0.25 for J 2 = |5l| 



and it decreases with J 2 as shown in Fig. 10. These deviations from field theory 
are probably due to higher order corrections in 1/S which have not been studied 
analytically so far. 

In the spiral phase, it is necessary to use a different NLSM which is known for 
5 = |5^, [53J. The field variable is now an SO (3) matrix R. The Lagrangian density 
is 

C = ^Tr(R-RP ) - ^Tr(R-R' Pl ), (8) 



where c = S(l + y)VT^f/y, g 2 = 2^(1 + y)/(l - yj/S with l/y = 4J 2 , and P 
and Pi are diagonal matrices with diagonal elements (1, 1, 2y(l —y)J {2y 2 — 2y + 1)) 
and (1,1,0) respectively. Note that there is no topological term. Hence there is 
no apparent difference between integer and half-odd-integer spin chains in the spiral 
phase. A one-loop renormalization group []52 and large N analysis |53| indicate that 
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the system should have a gap for all values of J2 and S, and that there is no reason 
for a particularly small gap at any special value of J2 . 

Finally, in the colinear phase, the NLSM is known for 8 — 1, i.e., for the spin 
ladder. The Lagrangian is the same as in (|7|), with c = 4Sy J 2^2 + 1), g 2 = 

J\ + 1/ J2/S and 9 = 0. There is no topological term for any value of S, and the 
model is therefore gapped. 

The field theories for general 5 in both the spiral and colinear phases are still 
not known. Although the results are qualitatively expected to be similar to the 
5 = case in the spiral phase and the 5 = 1 case in the colinear phase, quantitative 
features such as the dependence of the gap on the coupling strengths will require 
the explicit form of the field theory. 



4.2 Bosonization 

Another field theoretic method for studying spin systems in one dimension is the 
technique of bosonization E5L H7L BH ISBfl. This technique consists of mapping 



bosonic operators into fermionic ones, and then using whichever set of operators 
is easier to compute with. For instance, consider a model with a single species of 
fermion with a linear dispersion relation E(k) = ±vk, where the ± denotes the 
right- and left-moving fermions respectively (with the corresponding fields being 
denoted by ipR an d i/il), an d v denotes the velocity. Similarly, consider a model 
with a single species of boson with the dispersion relation E(k) = v\k\; the right- 
and left-moving fields are denoted by (f)R and 4>l respectively. Then it can be shown 
that these operators are related to each other as 

1 



4> 



-i2y^7T(f> R 



V2 



% L ~ -L= e t2 ^ L . (9) 
y/Zira 

The length parameter a is a cut-off which is required to ensure that the contri- 
bution from high-momentum modes do not produce divergences when computing 
correlation functions. It is convenient to define the two bosonic fields 

<j> = 4>r + 4> L , 

= - h + k- (10) 
Then the fermionic density is given by 

P - PO = ^R^R + $3h = ~ -j= , (11) 
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where p is the background density; fluctuations around this density are described 
by the quantum fields ip or <f>. 

Although the dispersion relation is generally not linear for all the modes of a 
given system, it often happens that the low-energy and long-wavelength modes can 
be studied using bosonization. For a fermionic system in one dimension, these 
modes are usually the ones lying close to the two Fermi points with momenta ±fcp 
respectively. One can define right- and left-moving fields ipn and ipL which vary 
slowly on the scale length a, 

${x) = $ R {x) e ikFX + i) L {x) e~ ikFX . (12) 

Quantities such as the density will generally contain terms which vary slowly as well 
as terms varying rapidly on the scale of a, 

P-Po = W = $ R j>R + Ml + e- m * x $ R j> L + e l2kFX ^ R 

One can compute various correlation functions in the bosonic language. Consider 
an operator of the form 

O m ,n = e l2 ^ im * +n§) . (14) 

We find the following result for the two-point equal-time correlation function at 
spatial separations which are much larger than the microscopic lattice spacing a, 

~2 is i ^2 ) 



(0| Td m!n (x)6l, ;n ,(0) |0> ~ 6 mm ,6 nn , (-) 2(m K+n ' K) , (15) 



where K denotes an interaction parameter which will be described below. Note that 
the correlation function decays as a power law. In the language of the renormaliza- 
tion group, the scaling dimension of O m , n is given by m 2 K + n 2 /K. 

We can now study a quantum spin chain using bosonization. To be specific, let 
us consider a spin- 1/2 chain described by the anisotropic Hamiltonian 

N J . . 

H = E [ o {S?S»i + S rS? + i) + JzS?S! +1 - hS? ] , (16) 

where the interactions are only between nearest neighbor spins, and J > 0. = 
Sf + iSf and S~ = S x — iSf are the spin raising and lowering operators, and 
h denotes a magnetic field. Note that the model has a (7(1) invariance, namely, 
rotations about the S z axis. When J z = J and h — 0, the U(l) invariance is 
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enhanced to an SU(2) invariance, because at this point the model can be written 
simply as H = JYsi&i' Si+i- Although the model in (|16j) can be exactly solved 
using the Bethe ansatz, and one has the explicit result that the model is gapless for 
a certain range of values of J z /J and h/J (see Ref. |59|), it is not easy to compute 
explicit correlation functions in that approach. We therefore use bosonization to 
study this model. 

We first use the Jordan- Wigner transformation to map the spin model to a model 
of spinless fermions. We map an f spin or a j spin at any site to the presence or 
absence of a fermion at that site. We introduce a fermion annihilation operator ipi 
at each site, and write the spin at the site as 

SI = $^-1/2 = 7^-1/2 

S~ = (-l)'Ae^«* , (17) 

where the sum runs from one boundary of the chain up to the (i — l) th site (we 
assume open boundary conditions here for convenience), = or 1 is the fermion 
occupation number at site i, and the expression for is obtained by taking the 
hermitian conjugate of , The string factor in the definition of S~ is added in order 
to ensure the correct statistics for different sites; the fermion operators at different 
sites anticommute, whereas the spin operators commute. 

We now find that 

H = - [ 77 ($A+i + h.c.) + J z (n i -l/2)(n i+1 -l/2) - h {hi -1/2) } . (18) 
i 1 

We see that the spin-flip operators Sf lead to hopping terms in the fermion Hamilto- 
nian, whereas the SfS* +l interaction term leads to an interaction between fermions 
on adjacent sites. 

Let us first consider the noninteracting case given by J z = 0. By Fourier trans- 
forming the fermions, tjjk = J2j ^je~' lk ^ a / yN, where a is the lattice spacing and the 
momentum k lies in the first Brillouin zone —it /a < k < it /a, we find that the 
Hamiltonian is given by 

H = "i< Wk , (19) 
k 

where 

u>k = — J cos(fca) — h . (20) 

The noninteracting ground state is the one in which all the single-particle states 
with ojk < are occupied, and all the states with ujj~ > are empty. If we set the 
magnetic field h — 0, the magnetization per site m = Y^iSf/N will be zero in the 
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ground state; equivalently, in the fermionic language, the ground state is precisely 
half-filled. Thus, for m — 0, the Fermi points (ojk = 0) lie at ka = ±7r/2 = k F a. Let 
us now add the magnetic field term. In the fermionic language, this is equivalent 
to adding a chemical potential term (which couples to hi or Sf). In that case, the 
ground state no longer has m = and the fermion model is no longer half-filled. 
The Fermi points are then given by ±fcp, where 

k F a = Ti (m + -) . (21) 

It turns out that this relation between kp (which governs the oscillations in the 
correlation functions as discussed below) and the magnetization m continues to 
hold even if we turn on the interaction J z , although the simple picture of the ground 
state (with states filled below some energy and empty above some energy) no longer 
holds in that case. 

In the linearized approximation, the modes near the two Fermi points have the 
velocities dujk/dk = ±v, where v is some function of J, J z and h. Next, we introduce 
the slowly varying fermionic fields ipR and ^ as indicated above; these are functions 
of a coordinate x which must be an integer multiple of a. Finally, we bosonize these 
fields. The spin fields can be written in terms of either the fermionic or the bosonic 
fields. For instance, S z is given by the fermion density as in Eq. ( pi) which then 
has a bosonized form given in Eq. ([TBI). Similarly, 

S + (x) = {-l) x ' a [ e - ikFX/a ft R (x) + e ikFX/a 4>l(x)] x 

[ e i7r dx\i,\x')^{x')+l/2a) + ^ c J ^ 22 ^ 

where ( — \) x / a = ±1 since x/a is an integer. This can now be written entirely in the 
bosonic language; the term in the exponential is given by 

f dx' ft(x')j>(x') = - 4= r dx ' in = ~7= [Mx)+h{x)\ , (23) 

J-oo y 7T J-oo ox 

where we have ignored the contribution from the lower limit at x' = — oo. 

We can now use these bosonic expressions to compute the two-spin, equal-time 
correlation functions G ab (x) =< 0\TS a (x)S b (0)\0 >. We find that 

G zz (x) 
G + -(x) + G-+(x) 



m 2 + -f- + c 2 



Ci cos(2k F x) 



(-l) x / a (-1)*/° cos(2A; F x) 

C 3 1/2K + ° 4 r 2K+l/2K ' ^ 
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where c±, C4 are some constants. The parameters K and v are functions of J z /J 
and h/J; the exact dependence can be found from the web site given in Ref. [ |39| . 
For h = 0, K is given by the analytical expression 

1 = 1 + -sin- 1 A. (25) 

K 71 J 

Note that at the SU(2) invariant point J z = J and h = 0, we have K = 1/2, and 
the two correlations G zz and G + ~ have the same forms. 

In addition to providing a convenient way of computing correlation functions, 
bosonization also allows us to study the effects of small perturbations. For instance, 
a physically important perturbation is a dimerizing term 

V = S £ (-1)* [ J - {S+Sr +1 + S'rS'+J + J Z S?S* +1 ] , (26) 

i 

where 5 is the strength of the perturbation. Upon bosonizing, we find that the 
scaling dimension of this term is K. Hence it is relevant if K < 2; in that case, it 
produces an energy gap in the system which scales with S as 

AE ~ d 1 '^ . (27) 

This kind of phenomenon occurs in spin-Peierls systems such as CuGeC>3; below a 



transition temperature T sp , they go into a dimerized phase which has a gap [57 



5 Density Matrix Renormalization Group Method 

The method which held the promise of overcoming the difficulty of exploding dimen- 
sionalities is the renormalization group (RG) technique, in which one systematically 
throws out the degrees of freedom of a many-body system. While this technique 
found dramatic success in the Kondo problem f58|| , its straightforward extension to 



interacting lattice models was quite inaccurate [59|. 

In early 1992, the key problems associated with the failure of the old RG method 
were identified and a different renormalization procedure based on the eigenvalues 
of the many-body density matrix of proper subsystems was developed [^, . This 
method has come to be known as the density matrix renormalization group (DMRG) 
method and has found dramatic success in solving quasi-one-dimensional many-body 
Hamiltonians. 

In a real-space RG approach, one begins by subdividing the total system into 
several blocks A n and proceeds to iteratively build effective blocks so that at each 
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iteration, each effective block represents two or more blocks of the previous iteration, 
without increasing the Fock space dimensionality of the blocks from what existed at 
the previous iteration. Usually, one starts with each A n consisting of a single site. 
Since the Hilbert space grows exponentially with the increase in system size, one 
truncates the number of states kept at each iteration. 

The main reason for the failure of the old RG methods is the choice of the states 
retained at each stage of the iteration |60|| . White [|6T]| , recognized that the weakness 



of the old RG procedure was in the truncation of the Fock space of a block based 
on the eigenvalues of the block Hamiltonian being renormalized. He replaced this 
choice by introducing a completely different truncation scheme than the one that 
was used in the old quantum RG procedures. The choice is the eigenvalues of the 
reduced density matrix of the block constructed from the desired state of the full 
Hamiltonian. The truncated Fock space is now spanned by the m eigenvectors of 
the reduced density matrix of order I x I (m < I) corresponding to the m highest 
eigenvalues of the reduced many-body density matrix. The reason for choosing the 
eigenvalues of the reduced density matrix as a criterion for implementing a cut- 
off is that, the larger the density matrix eigenvalue, the larger is the weight of 
the eigenstate of the density matrix in the expectation value of any property of the 
system. This result becomes evident when all the dynamical operators are expressed 
as matrices in the basis of the eigenvectors of the density matrix. The expectation 
value of any operator A is simply 

< A > = J2 A HPi I Eft > ( 28 ) 

i i 

where pi is the density matrix eigenvalue. The larger the value of a particular pi, 
the larger is its contribution to the expectation value, for a physically reasonable 
spread in the diagonal matrix elements An. 

The many-body density matrix of a part of the system can be easily constructed 
as follows. Let us begin with given state \tp >s of S, which is called the universe or 
superblock, consisting of the system (which we call a block) A and its environment 
A'. Let us assume that the Fock space of A and A' are known, and can be labeled 
as \i >a and \j >a> respectively. The representation of \ip >$ in the product basis 
of %a and %a' can be written as 

l^>5 = Y f i>i&>A*\i>A! > ( 29 ) 

ij 

where we assume the coefficients ipij to be real, without loss of generality. Then the 
reduced many-body density matrix for block A, is defined as 

Pki = ^24>kjipij ■ (30) 
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The eigenvalue pi of the density matrix p gives the probability of finding the cor- 
responding eigenstate \pi >a in the projection of \tp >s on block A. It therefore 
follows that the eigenvectors with highest eigenvalues of the density matrix of A, are 
the optimal or most probable states to be retained while the system is augmented. 

In the early literature in quantum chemistry, the eigenvectors corresponding to 
large eigenvalues of one-particle density matrices were employed as the orbital basis 
for carrying out a configuration interaction (CI) calculation. The eigenvectors of the 
density matrix were called the 'natural' orbitals, and it was observed that the CI 
procedure converged rapidly when the 'natural' orbitals were employed in setting 
up the Slater determinants H62]] . 

The DMRG scheme differs from the 'natural' orbital scheme in two important 
respects: (i) the reduced density matrices are many-body density matrices, and (ii) 
the size of the system in terms of the number of sites being studied at each itera- 
tion is usually augmented by two sites. However, the Hamiltonian matrix that one 
encounters from iteration to iteration, remains roughly of the same order while the 
matrix elements keep changing (renormalized) . In this sense, the procedure can be 
called a renormalization procedure. The coupling constants (the Hamiltonian ma- 
trix elements) keep changing while the system size increases, as in the RG procedure 
carried out within a blocking technique. 

5.1 Implementation of the DMRG method 

We now describe the procedure to carry out the computations. One starts the 
computation with a small size system, 2n, which can be exactly solved, 1 < n < 4, 
depending on the degree of freedom at each site. By exact diagonalization, one 
gets the desired eigenstate of that system. The density matrices of the left and 
right blocks, each consisting of n sites (in principle it is not necessary to have the 
same number of sites for the two blocks, although in practice this is what is most 
generally used), are obtained from the desired eigenstate. The density matrices 
are diagonalized and at the first iteration usually all the density matrix eigenvectors 
(DMEV) are retained. The Hamiltonian matrix of the left and right blocks (denoted 
by A and A') obtained in any convenient basis are transformed into the density 
matrix eigenvector basis. So also are the matrices corresponding to the relevant site 
operators in both blocks. Now, the iterative procedure proceeds as follows. 

1. Construct a superblock S = A • »A', consisting of the block A, two additional 
sites a and o' and the block A'. Thus, at the first iteration, the system S has 
n+1 + 1 +n — In + 2 sites. 
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2. Set up the matrices for the total Hamiltonian of the superblock S in the 
direct product basis of the DMEV of the blocks A and A' and the Fock space 
states of the new sites. Considering that the new sites are spin-S 1 sites with 
(2S + 1) Fock states each, the order of the total Hamiltonian matrix will be 
m 2 (2S + l) 2 x m 2 (2S + l) 2 , where m is the dimension of the block DMEV 
basis. 

3. Diagonalize the Hamiltonian of the superblock S = 2n + 2 to find the desired 
eigenstate \ip >. Using the state \ip >, evaluate all properties of the superblock 
of interest. 

4. Construct the reduced many-body density matrix, p, for the new block A*. If 
the system does not possess reflection symmetry, construct the density matrix, 
p', for the new right block mA' as well. 

5. Diagonalize the density matrix, p, and if necessary p'. Usually, the density 
matrix is block-diagonal in the z-component of the total spin of the block, 
and it becomes computationally efficient to exploit such quantum numbers. 
Construct a nonsquare matrix O, with m columns, each column being an 
eigenvector of the density matrix corresponding to one of the m largest eigen- 
values. The number of rows in the matrix O corresponds to the order of the 
density matrix. 

6. Construct the matrices corresponding to the Hamiltonian, Ha», of the new left 
block A», and the site spin operators (S z and S + ) of all the necessary sites. 
The S~ operators are simply the adjoints of the S + operators. 

7. Renormalize all the matrices corresponding to the block and site operators by 
using the RG transformation matrix O, e. g., Ha, = OHj^O* . The resulting 
renormalized matrices are of order m x m and the procedure amounts to a 
simultaneous change of basis and a truncation. 

8. Replace the block A by A*. If the system does not possess reflection symmetry, 
replace A' by 

9. Go to step 1. 

Using the block-diagonal nature of the density matrix, besides reducing the re- 
quirement in CPU time, also allows one to label the DMEV by the appropriate 
^-component of the total spin of the block (M SjA ). The Fock space of the individual 
sites that are added at each iteration are eigenstates of the site spin and number 
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operators. This allows us to target a definite projected spin (M s ) state of the total 
system. 

We now briefly describe the mathematical notations that we have used so far for 
various states. A state of A* is given by the tensor product of a state of A with 
quantum number q and an index i, and a state o of the additional site. Thus, 

\q,i,a> A . = \q,i> A x\cr> . (31) 

A state of a superblock 5* = A • »A' is given by 

\q A ,li,a;qA',v,T > = \q A , fi,a > A , x\q A ,,u,r >, A > . (32) 

The eigenstate of the Hamiltonian of the super-block can be written as 

l^>5 = E ^ qA " a ' T \qA^,cr;q A ,,v,T> s . (33) 

The density matrix for A* then will have a block structure and can be expressed as 

pir = e €T" C,T ^ A " a ' T ■ (34) 

The above algorithm is called the infinite lattice DMRG algorithm because this 
procedure is best suited for the system in the thermodynamic limit, i.e., when the 
properties of the system are extrapolated to the infinite system size limit. 



5.2 Finite size DMRG algorithm 

If we are interested in accurate properties of the system at a required size, then 
it is possible to improve upon the accuracies obtainable from the infinite DMRG 
procedure. This involves recognizing that the reduced many-body density matrices 
at each iteration correspond to a different system size. For example, when we carry 
out the DMRG procedure to obtain the properties of a system of 2M sites, at an 
iteration corresponding to 2p sites (n < p < M), the reduced density matrix we 
construct is that of a block of p sites in a system of 2p sites. However, if our interest 
is in the 2M-site system, we should employ the density matrix of the block of p 
sites in a 2M-site system. It is possible to construct, iteratively, the p-site reduced 
density matrix of the 2M-site system. This is achieved by the so called finite-size 



algorithm [5T . This method provides highly accurate solutions even when the states 



of the full Hamiltonian have inhomogeneous (symmetry breaking) properties. 
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To obtain the 2M-site result, we should perform the infinite lattice algorithm up 
to I = (M — 1) sites first storing all operators in each iteration. Now the algorithm 
for finite lattices with reflection symmetry (left block = right block), proceeds as 
follows. 

1. On reaching a system size of 2M sites, obtain the density matrix of the block 
of M sites. 

2. Use the density matrix of M sites on the left and that of (M — 2) sites on the 
right, add two new sites as in the infinite DMRG procedure and obtain the 
desired eigenstate of the 2M system. 

3. Now obtain the reduced density matrix of the (M + l) sites from the eigenstate 
of the previous iteration obtained in the direct product basis of the DMEVs 
of the M-site, (M — 2)-site density matrices, and the Fock space states of the 
individual sites. 

4. Go back to step 2, replacing M by (M + l) and (M — 2) by (M — 3) and iterate 
until a single site results on the right and (2M — 3) sites result on the left. 

5. Since the system has reflection symmetry, use the density matrix of the [2M — 
3) sites on the right and construct the 2M system as built-up from three 
individual sites on the left and (2M — 3) sites on the right. Obtain the desired 
eigenstate of the 2M system in this basis. 

6. Now obtain the new 2-site density matrix on the left and (2M — 4) site density 
matrix on the right. Replace the single-site on the left by two sites and (2M— 3) 
sites on the right by (2M — 4) sites in step 5. 

7. Repeat steps 5 and 6 until (M — 1) sites are obtained both on the left and 
right. The properties of the 2M system obtained from the eigenstates at this 
stage corresponds to the first iteration of the finite-size algorithm. We can 
now go back to step 1 and carry through the steps to obtain properties at 
later iterations of the finite-size DMRG algorithm. 

In systems without reflection symmetry, the DMEVs of the right and left parts 
are not identical even if the sizes of the reduced systems are the same. The finite- 
DMRG algorithm in this case involves first constructing the density matrices of the 
left part for sizes greater than M and on reaching the density matrix of (2M — 3) 
sites, reducing the size of the left-part and increasing that of the right, from one site 
to (M — 1). This will result in the refined density matrices of both the right and 
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the left block of the total system, for block sizes of (M — 1). At this stage, we can 
compute all the properties and continue the reverse sweep until the right block is 
of size (2M — 3) and the left block is of size 1. The forward sweep that follows will 
increase the block size on the left and decrease that of the right. We would have 
completed the second iteration when the two block sizes are equal. The forward 
and reverse sweeps can be continued until we reach the desired convergence in the 
properties of the whole system. 



5.3 Calculation of properties in the DMRG basis 

At the end of each iteration, one can calculate the properties of the targeted state 
|63| . The reduced many-body density matrix computed at each iteration can be 
used to calculate the static expectation values of any site operator or their products. 
Care should be taken to use the density matrices appropriate to the iteration. The 
expectation value of a site property corresponding to the operator A, can be written 
as 

< A > = Tr (pAi) . (35) 

p is the density matrix of the block in which the site i is situated and Aj is the matrix 
of the renormalized site operator at site i. For calculating correlation functions, one 
can use a similar equation. The correlation function between two site operators 
belonging to separate blocks can be written as 

< AiAj > = Tr(pAjAj) . (36) 

However, the accuracy of this procedure turns out be very poor if the sites i and 
j belong to the same block ||61|| . The reason is that a feature implicit in the above 
procedure is the resolution of identity by expansion in terms of the complete basis. 
Unfortunately, the basis in which the site operators are represented is incomplete 
and such an expansion is therefore error prone. To circumvent this difficulty, it 



has been suggested |61[ that, one obtains the matrix representation of the products 
of the site operators from the first occurrence of the product pair < ij > and, by 
renormalizing the product operator A^Aj, at every subsequent iteration until the 
end of the RG procedure. Then, the correlation function between Aj and Aj (where 
i and j belong to the same block) can be evaluated as 

< A^ > = Tr(pAiAj) . (37) 

This procedure is found to be more accurate in most cases. 
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5.4 Remarks on the applications of DMRG 



The DMRG method is currently the most accurate method for large quantum lattice 
models in one dimension. It can be applied to interacting bosonic, fermionic or spin 
models as well as to models which have interactions amongst them. The overall 
accuracy of the DMRG method is exceptionally high for one-dimensional systems 
with only nearest neighbor interactions. For a spin- 1/2 chain where exact Bethe 
ansatz ground-state energy is available, the DMRG ground state energy per site in 
units of the exchange constant J, is found to be accurate to seven decimal places 



with a cut-off m = 100 ||61|| . The method is found to be almost as accurate for the 
one-dimensional Hubbard model, where again it is possible to compare the DMRG 
results with exact results obtained from the exact analytic Bethe ansatz solution 

0- 



Since higher dimensionality is equivalent to longer-range interactions within one 
dimension, the model also restricts the range of interactions in one dimension. It has 
been noted that the number of DMEVs that should be retained in a calculation on 
higher dimensional systems, for accuracies comparable to accuracy in one dimension, 
scales exponentially with dimensionality. Thus, to obtain accuracy comparable to 
that obtained in a chain of L-sites for a cut-off m, in a L x L square lattice, the 
number of DMEVs needed to be retained for the corresponding 2-dimensional lattice 
is ~ m 2 . 



Extending the range of interactions to next-nearest neighbors does not signifi- 



cantly deteriorate the accuracy J65|. However, inclusion of cyclic boundary condi- 
tions reduces the accuracy of the method significantly, although in one- dimension, 
the DMRG method still would outperform any other method for the same system 
size. In the DMRG procedure, the most accurate quantity computed is the total 
energy. In dealing with other quantities such as correlation functions, caution must 
be exercised in interpreting the results. 

The density matrix eigenvalues sum to unity and the truncation error, which is 
defined as the sum of the density matrix eigenvalues corresponding to the discarded 
DMEV, gives a qualitative estimate as to the accuracy of the calculation as well as 
providing a framework for the extrapolation to the m — > oo limit. The accuracy of 
the results obtained in this way is unprecedented |66], ^7] . The accuracy of the ground 



state energy per site for the spin-1 chain is limited by the precession of machine 
arithmetic, i.e., e$ = 1.401484038971(4). Similarly, the accuracy persists even while 
calculating for the Haldane gap, e.g., the gap is evaluated to be 0.41050(2). 

Another aspect of the DMRG technique worth noting is that the method is 
best suited for targeting one eigenstate at a time. However, it is possible to obtain 
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reasonable results for a set of states by using an average many-body reduced density 
matrix constructed as a weighted sum of the density matrices corresponding to each 
of the states in question. One way of constructing the average density matrix is by 
using a statistical weight for the chosen set of states; the averaged density matrix 
in this instance is given by 

Pp-M = ^Mkj Vtyy exp[-/3ej] / exp[-/3ei] , (38) 

ij i 

where (3 = 1/ksT with ks and T being the Boltzmann constant and temperature 
respectively. One can thus extend the DMRG method to finite temperatures. 

Finite size algorithms have been used extensively to study the edge states and 
systems with impurities, where substantial improvement of the accuracy is needed 
to characterize the various properties of a finite system. The DMRG method has 
been applied to diverse problems in magnetism: study of spin chains with s > 1/2 
chains with dimerization and/or frustration p], BR p9l [70], coupled spin chains 



[65 , 71, 72 1, to list a few. Highly accurate studies of the structure factor and string 



order parameter (topological long range order) [EJJ as well as edge states in Haldane 
phase systems ]FJJ have been performed. Dynamical properties for both spin and 
fermionic systems with DMRG have also been reported within the Maximum En- 
tropy Method J74| as well as the continued fraction |75[] and correction vector [76| ap- 



proaches. Also, DMRG has been successfully formulated to obtain low-temperature 
thermodynamic properties for various spin systems [[T_7], [TJJ, and the solution of mod- 
els of spin chains dynamically coupled to dispersionless phonons [[Fj|. Additionally, 



Nishino and Okunishi have derived two reformulations of DMRG, namely, the prod- 
uct wave function renormalization group (PWFRG) []80j, and the corner transfer 
matrix renormalization group (CTMRG) methods. These methods offer the 
means of calculating dynamical correlation functions in spin chains as well as highly 
accurate results for the 2-dimensional Ising model at criticality. 



6 Frustrated and Dimerized Spin Chains 

It is well known that the one-dimensional XY chain can be mapped on to a one- 
dimensional noninteracting spinless fermion model. The isotropic spin chain will 
then map on to a chain of interacting spinless fermions. According to Peierls' theo- 
rem, a partly filled one-dimensional band of noninteracting fermions is unstable with 
respect to a lattice distortion that results in an insulating ground state. It has been 
shown that introduction of interactions in the Peierls system leads to a stronger 
instability. The mapping between the Heisenberg spin chains with equal nearest 
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neighbor exchange interactions (uniform spin chain) and the spinless fermion model 
suggests that such a spin chain is also unstable with respect to a lattice distortion 
leading to alternately strong and weak nearest neighbor exchange constants, i.e., a 
dimerized spin chain. What is of importance is the fact that such a dimerization is 
unconditional: no matter how strong the lattice is, the lattice dimerizes, since the 
exchange energy gained due to dimerization always exceeds the strain energy. This 
result follows from the fact that the gain in exchange energy varies as <5 2 ln<5 while 
the strain energy loss varies as 5 2 , where 5 is the magnitude of dimerization that 
leads to the nearest neighbor exchange constants alternating as J(l±5). 

In recent years, many systems which closely approximate the one-dimensional 
spin chain have been synthesized. What has been observed in these experimental 
systems is that besides the nearest neighbor antiferromagnetic exchange, there also 
exists a second neighbor exchange J2 of the same sign and comparable magnitude. 
Such a second neighbor interaction has the effect of frustrating the spin alignment 
favored by the nearest neighbor interaction. Therefore, a realistic study of these 
systems require modeling them using both dimerization and frustration. Theoreti- 
cally, spin chains with only frustration ( J — J 2 model) were studied by Majumdar 
and Ghosh. Interestingly, they showed that for J 2 = J/2, the ground state is doubly 
degenerate and is spanned by the two possible Kekule structures (Fig. 11). It is 
quite gratifying to note that a century after the Kekule structure for benzene was 
proposed, there actually exists a Hamiltonian for which Kekule structure happens 
to be the ground state! 

While most of the discussion above is restricted to spin-1/2 chains, there has been 
considerable interest in the higher-spin chains following the conjecture of Haldane 
which predicts that for uniform spin chains, the excitation spectrum of integer spin 
chains is qualitatively different from that of half-odd-integer spin chains. The latter 
have a gapless excitation spectrum while the excitation spectrum of the former are 
gapped. The synthesis and study of integer spin chains have indeed confirmed this 
conjecture. 

Notwithstanding many interesting exact analytic solutions for spin chains, there 
still exist a large number of situations for which such solutions have been elusive. 
The exact solutions are basically confined to the uniform Heisenberg model and 
the frustrated and dimerized model along the line 2J 2 + 8 = 1 in the J 2 — 5 
plane, with J = 1. Reliable numerical study of these models therefore requires 
developing techniques which are highly accurate so that the results of large finite 
systems can be scaled or extrapolated to the thermodynamic limit. As has already 
been discussed, the DMRG technique is ideally suited due to its high accuracy for 
quasi-one-dimensional systems. 
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The Hamiltonian for the frustrated and dimerized spin chain is given in Eq. 
and is schematically shown in Fig. 8. A few low-lying states in a sector with a given 
value of the total spin component, Ms are obtained in representative points in the 
J 2 — 5 plane, using the DMRG method. The ground state is always the first (lowest 
energy) state in the Ms = sector. The accuracy of the DMRG method depends 
crucially on the number of eigenstates of the density matrix, m, which are retained. 
Working with m = 100 to 120 over the entire J 2 — 5 plane gives accurate results. 
This can be verified by comparing the DMRG results for these m values with exact 
numerical diagonalizations of chains with up to 16 sites for spin-1 systems |82| and 
22 sites for spin-1/2 systems [33|j. The chain lengths studied vary from 150 sites for 
J 2 > to 200 sites for J 2 = 0. The DMRG results are also tracked as a function of 
N, the chain length, to verify that convergence is reached well before 150 sites in all 
cases. The numerical results are much better convergent for open chains than for 
periodic chains, a feature generic to the DMRG technique [^T| ^7] . 



The quantum phase diagrams obtained for a spin-1/2 chain is shown in Fig. 12. 
The system is gapless on the line A running from J 2 = to J 2c = 0.241 for 5 = 0, 
and is gapped everywhere else in the J 2 — 5 plane. There is a disorder line B given 
by 2J 2 + 5 = 1; the peak in the structure factor S(q) is at q max = ir to the left of 
B (region I), decreases from tt to 7r/2 as we go from B up to the line C (region II), 
and is at q m ax — 7r /2 to the right of C (region III). This is in agreement with the 
results obtained in section 4 using the NLSM approach. The correlation length £ 
goes through a minimum on line B. 

In the spin-1 case (Fig. 10), the phase diagram is more complex. There is a solid 
line marked A which runs from (0, 0.25) to about (0.22 ±0.02, 0.20 ±0.02) shown by 
a cross. To within numerical accuracy, the gap is zero on this line and the correlation 
length £ is as large as the system size N. The rest of the 'phase' diagram is gapped. 
However the gapped portion can be divided into different regions characterized by 
other interesting features. On the dotted lines marked B, the gap is finite. Although 
£ goes through a maximum when we cross B in going from region II to region I or 
from region III to region IV, its value is much smaller than N. There is a dashed 
line C extending from (0.65, 0.05) to about (0.73, 0) on which the gap appears to be 
zero (to numerical accuracy), and £ is very large but not as large as N. The straight 
line D satisfying 2 J 2 + 5 = 1 extends from (0, 1) to about (0.432, 0.136). Regions II 
and III are separated by line E which goes down to about (0.39,0). Across D and 
E, the peak in the structure factor decreases from it (Neel) in regions I and II to less 
than tt (spiral) in regions III and IV. In regions II and III, the ground state for an 
open chain has a four- fold degeneracy (consisting of states with £7 = and £7 = 1), 
whereas it is nondegenerate in regions I and IV with £7 = 0. The regions II and III, 
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where the ground state is four-fold degenerate for an open chain, can be identified 
with the Haldane phase; the regions I and IV correspond to a non-Haldane singlet 
phase. The lines B, D and E meet in a small region V where the ground state of 
the system is numerically very difficult to find. Note that the numerically zero gap 
at (0.73,0) is unexpected from either bosonic mean-field theory ]84| or the NLSM 
approach discussed earlier. 

For the spin-1 system, there is a striking similarity between the ground state 
properties of the dimerized and frustrated model (f|) as a function of J 2 (with 5 = 0) 
and the biquadratic model, 

H = Y, [ S* • Sm + /? ( Si ■ S m f ] , (39) 



as a function of (positive) (3 |59|. For J 2 < 0.39 and (3 < 1/3, both models are in the 



Neel phase and are gapped. For J 2 > 0.39 and /5 > 1/3, the two models are in the 
spiral phase and are generally gapped; however, model (f|) is 'gapless' for J 2 = 0.73 
while model fl39|) is gapless for (3=1. Qualitatively, the cross-over from the Neel 
to the spiral phase (but not the gaplessness at a particular value of J 2 or (3) can be 
understood through the following classical argument. Let us set the magnitudes of 
the spins equal to 1 and define the angle between spins Sj and S i+n to be n8. The 
angle 8 can be obtained by minimizing cos 8 + J 2 cos 28 in (TJ), and cos# + (3 cos 2 8 
in (|39D . This gives us a Neel phase {8 = n) if J 2 < 1/4 and (3 < 1/2 in the two 
models, and a spiral phase for larger values of J 2 and (3 with 8 = cos _1 (— 1/4 J 2 ) 
and 8 = cos _1 (— 1/2(3) respectively. The actual crossover points from Neel to spiral 
is different for spin-1 than these classical values. In the classical limit S — > 00, the 
ground state of the model is in the Neel phase for 4J 2 < 1 — 5 2 , in a spiral phase for 
1 - 5 2 < 4J 2 < (1 - 5 2 )/5 and in the colinear phase for (1 - 5 2 )/5 < 4J 2 (Fig. 9). 

As can be seen from Fig. 8, 5 = 1 results in two coupled spin chains wherein 
the interchain coupling is 2 and the intrachain coupling is J 2 . By using DMRG, 
one can study the dependence of the gap A and the two-spin correlation function 
C(r) on the interchain coupling J. In Fig. 13 is plotted A vs J for both spin-1/2 
and spin-1 systems. For spin-1/2, the system is gapped for any nonzero value of the 
interchain coupling J, although the gap vanishes as J — > 0. The gap increases and 
correspondingly the correlation length decreases with increasing J. In the case of 
coupled spin-1 chains, one finds the somewhat surprising result that both the gap 
and the correlation length £ are fairly large for moderate values of J. Note that the 
variation of the gap with J for spin-1 (shown as circles) is much less than that for 
spin-1/2 (crosses). 

The NLSMs derived in section 4 can be expected to be accurate only for large 
values of the spin S. It is interesting to note that the numerically obtained 'phase' 
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boundary between the Neel and spiral phases for spin-1 is closer to the the classical 
(5 — > oo) boundary 4J 2 = 1 —S 2 than for spin-1/2. For instance, the crossover from 
Neel to spiral occurs, for 8 = 0, at J 2 = 0.5 for spin-1/2, at 0.39 for spin-1, and at 
0.25 classically. 

To conclude this section, we have studied a two-parameter 'phase' diagram for 
the ground state of isotropic antiferromagnetic spin-1/2 and spin-1 chains. The spin- 
1 diagram is considerably more complex than the corresponding spin-1/2 chain, with 
surprising features like a 'gapless' point inside the spiral 'phase'; this point could 
be close to a critical point discussed earlier in the literature |5U], |55|| . It would be 
interesting to establish this more definitively. Our results show that frustrated spin 
chains with small values of S exhibit features not anticipated from large S field 
theories. 



7 Alternating (Si, 62) ferrimagnetic spin chains 

Ferrimagnets belong to a class of magnets which show spontaneous magnetization 
below a certain critical temperature. There have been a number of experimen- 
tal efforts to synthesize molecular materials showing spontaneous magnetization 
at low temperatures [^6], 87]. These are quasi-one-dimensional bimetallic molecu- 
lar magnets in which each unit cell contains two spins with different spin values, 
with the general formula ACu(pbaOH)(H 2 0)3.2H 2 0, where pbaOH is 2-hydroxo- 
l,3-propylenebis(oxamato) and A = Mn, Fe, Co, Ni; they belong to the alternating 
or mixed spin chain family |57|, |8l| . These alternating spin compounds have been 



seen to exhibit ferrimagnetic behavior. It has been possible to vary the spin at 
each site from low values where quantum effects dominate to large values which are 
almost classical. 

The thermodynamic behavior of these alternating spin compounds is very in- 
teresting [RR R!J. In very low magnetic fields, these systems show one-dimensional 
ferrimagnetic behavior. The \T vs T (where \ is the magnetic susceptibility and T 
the temperature) plots show a rounded minimum. As the temperature is increased, 
XT decreases sharply, goes through a minimum before increasing gradually. The 
temperature at which this minimum occurs differs from system to system and de- 
pends on the site spins of the chain. The variation of the field induced magnetization 
with temperature is also interesting as the ground state is a magnetic state. These 
exciting observations have motivated us to study the ferrimagnetic systems with 
arbitrary spins s\ and s 2 alternating from site to site. It would also be of interest 
to know the thermodynamic properties of these systems with varying s± and s 2 . 
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7.1 Ground state and excitation spectrum 

We start our discussion with the Hamiltonian for a chain with spins Si and s 2 on 
alternating sites (with si > s 2 , without loss of generality). 

H = J [ C 1 + 5 ) Si,n ■ S 2 ,„ + (1-5) S 2 , n ■ S ljn+1 ] , (40) 

n 

where the total number of sites is 2N and the sum is over the total number of unit 
cells N. Si tn corresponds to the spin operator for the site spin in the n-th unit 
cell. The exchange integral J is taken to be positive for all our calculations; S is the 
dimerization parameter and it lies in the range {0, 1}. 

Before describing our numerical results, we briefly summarize the results of a 
spin wave analysis for the purposes of comparison f78fl . We will first state the 
results for 5 = 0. According to spin wave theory, the ground state has total spin 
Sq = N(si — s 2 ). Let us define a function 

oo{k) = J ^{ Sl -s 2 ) 2 + 4 Sl s 2 sm 2 (k/2) , (41) 

where k denotes the wave number. Then the ground state energy per site is given 
by 

En 1 r dk , , , . , . . . . 

eo = ^ = - Js lS2 + - — - J(si + s 2 ) + u(k) ] . (42) 

ZiV Z JO 7T 

The lowest branch of excitations is to states with spin S = Sq—1, with the dispersion 

wi{k) = J(-s 1 + s 2 ) + to(k) ; (43) 

the gap vanishes at k = 0. There is a gapped branch of excitations to states with 
spin S = Sq + 1, with the dispersion 

u 2 (k) = J (si — s 2 ) + u{k) ■ (44) 

the minimum gap occurs at k = and is given by A = 2J(s\ — s 2 ). In the ground 
state with S z = So, the sublattice magnetizations are given by the expectation 
values 

< Sl,n > 

< $2,71 > 

The various two-spin correlation functions decay exponentially with distance; the 
inverse correlation length is given by £ -1 = ln(si/s 2 ). The results with dimerization 
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(5 > 0) are very similar. In fact, within spin wave theory, the minimum gap A to 
states with spin S = Sq + 1 is independent of 5. 

We now use the powerful DMRG method to study the system defined by Eq. ((40|) 
both with and without dimerization, 5^0 and 5 = respectively. We have consid- 
ered alternating spin-3/2/spin-l (hereafter designated as (3/2, 1)), spin-3/2 /spin- 1/2 
(to be called (3/2, 1/2)) and spin-1, 1/2 (called (1, 1/2)) chains with open boundary 
condition for the Hamiltonian ([H]). We compute the ground state properties for 
these three systems by studying chains with 80 to 100 sites. The number of domi- 
nant density matrix eigenstates, m, that we have retained at each DMRG iteration 
also varies between 80 to 100. With the increase of the Fock space dimensionality 
of the site spins, we increase m to obtain more accurate results. We follow the 
usual steps for the "infinite system" DMRG method discussed above |65|, |90 



except that the alternating chains studied here are not symmetric between the left 
and right halves; hence the density matrices for these two halves have to be sepa- 
rately constructed at every iteration of the calculations. We have also verified the 
convergence of our results by varying the values of m and the system size. 

The ground states of all the systems lie in the S z = N(si — s 2 ) sector, as 
verified from extensive checks carried out by obtaining the low-energy eigenstates 
in different S z sectors of a 20-site chain. A state corresponding to the lowest energy 
in S z = N(si — s 2 ) is found in all subspaces with \S Z \ < N(si — S2), and is absent 
in subspaces with IS' 2 ! > N(si — S2). This shows that the spin in the ground state 
is Sq = N(si — S2). (Actually, the lowest energy states in the different S z sectors 
are found to be degenerate only up to 10~ 5 J. Such small errors are negligible for 
studying thermodynamics at temperatures larger than, say, 10~ 2 J). 

In Fig. 14, we show the expectation value of site spin operator S z n ( spin density) 
at all the sites for the (3/2, 1), (3/2, 1/2) and (1, 1/2) chains. The spin densities are 
uniform on each of the sublattices in the chain for all the three systems. For the 
(3/2, 1) chain, the spin density at a spin-3/2 is 1.14427 (the classical value is 3/2), 
while, at a spin-1 site it is —0.64427 (classical value 1). For the (3/2, 1/2) case, the 
spin density at a spin-3/2 site is 1.35742 and at a spin- 1/2 site it is —0.35742. For 
the (1,1/2) case, the value at a spin-1 site is 0.79248 and at a spin-1/2 site it is 
—0.29248. These can be compared with the spin wave values of 1.040 and —0.540; 
1.314 and —0.314; and 0.695 and —0.195 for the spin-si and spin-s 2 sites of the 
(3/2,1); (3/2,1/2); and (1,1/2) systems respectively. We note that the spin wave 
analysis overestimates the quantum fluctuations in case of systems with small site 
spin values. We also notice that there is a greater quantum fluctuation when the 
difference in site spin \s\ — S2I is larger. This is also seen in spin- wave theory. The 
spin density distribution in an alternating (s!,s 2 ) chain is more similar to that of 
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a ferromagnetic chain rather than an antiferro magnet, with the net spin of each 
unit cell perfectly aligned (but with small quantum fluctuations on the individual 
sublattices). In a ferromagnetic ground state, the spin density at each site has 
the classical value appropriate to the site spin, whereas for an antiferro magnet, 
this averages out to zero at each site as the ground state is nonmagnetic. From 
this viewpoint, the ferrimagnet is similar to a ferromagnet and is quite unlike an 
antiferromagnet. The spin wave analysis also yields the same physical picture. 

Because of the alternation of spin-si and spin-s 2 sites along the chain, one has to 
distinguish between three different types of pair correlations, namely, < Si S( n >, 
< S^fiS^n > and < Si S% n >. We calculate all the three correlation functions with 
the mean values subtracted out, since the mean values are nonzero in all these three 
systems unlike in pure antiferromagnetic spin chains. In the DMRG procedure, 
we have computed these correlation functions from the sites inserted at the last 
iteration, to minimize numerical errors. In Fig. 15, we plot the two-spin correlation 
functions in the ground state as a function of the distance between the spins for 
an open chain of 100 sites for all three cases. All three correlation functions decay 
rapidly with distance for each of the three systems. From the figure it is clear that, 
except for the < Sl^S^ n > correlation, all other correlations are almost zero even 

for the shortest possible distances. The < Sl^S^n > correlation has an appreciable 
value [-0.2 for (3/2,1), -0.07 for (3/2,1/2) and -0.094 for (1,1/2)] only for the 
nearest neighbors. This rapid decay of the correlation functions makes it difficult 
to find the exact correlation length £ for a lattice model, although it is clear that 
£ is very small (less than one unit cell) for the (3/2, 1/2) and (1, 1/2) cases, and a 
little greater (1 < £ < 2) for the (3/2, 1) system. Spin wave theory gives £ = 2.47 
for (3/2,1), £ = 0.91 for (3/2,1/2), and £ = 1.44 for (1,1/2) cases. (We should 
remark here that our £ is not to be confused with the conventional definition of 
the correlation length; the latter is actually infinite in these systems due to the 
long-range ferrimagnetic order). 

The lowest spin excitation of all the three chains is to a state with S — Sq — 1. 
To study this state, we target the 2 nd state in the S z = S G — 1 sector of the chain. 
To confirm that this state is a S = Sg — 1 state, we have computed the 2 nd state in 
S z = sector and find that it also has the same energy. However, the corresponding 
state is absent in S z sectors with \S Z \ > Sg — 1. Besides, from exact diagonalization 
of all the states of all the Si — s 2 alternating spin chains with 8 sites, we find that 
the energy orderings of the states is such that the lowest excitation is to a state with 
spin S — Sq — 1- We have obtained the excitation gaps for all the three alternating 
spin chains in the limit of infinite chain length by extrapolating from the plot of 
spin gap vs the inverse of the chain length (Fig. 16). We find that this excitation is 
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gapless in the infinite chain limit for all three cases. 

To characterize the lowest spin excitations completely, we also have computed 
the energy of the S = Sq + 1 state by targeting the lowest state in the S z = So + 1 
sector. In Fig. 17, we plot the excitation gaps to the S = Sq + 1 state from 
the ground state for all three systems as a function of the inverse of the chain 
length. The gap saturates to a finite value of (1.0221 ± 0.0001) J for the (3/2, 1) 
case, (1.8558±0.0001)Jfor (3/2, 1/2), and (1.2795±0.0001) J for (1, 1/2). It appears 
that the gap is also higher when the difference in site spins \s\ — S2I, is larger. The 
site spin densities expectation values computed in this state for all three cases are 
found to be uniform (i.e., independent of the site) on each of the sublattices. This 
leads us to believe that this excitation cannot be characterized as the states of a 
magnon confined in a box, as has been observed for a spin-1 chain in the Haldane 



phase |56 



We have also studied the spin excitations in the dimerized alternating (s l5 S2) 
chains, defined in Eq. ((40|). We calculate the lowest spin excitation to the S = Sq — 1 
state from the ground state. We find that the S = Sg — 1 state is gapless from the 
ground state for all values of 5. This result agrees with the spin wave analysis of 
the general (si,S2) chain. The systems remain gapless even while dimerized unlike 
the pure antiferromagnetic dimerized spin chains. There is a smooth increase of the 
spin excitation gap from the ground state to the S = So + 1 state with increasing 
5 for all three systems studied here. We have plotted this gap vs 5 in Fig. 18. The 
gap shows almost a linear behavior as a function of S, with an exponent of 1.0 ±0.01 
for all three systems. This seems to be an interesting feature of all ferrimagnets. 
The spin wave analysis however shows that this excitation gap is independent of S 
for the general (si, sz) chain. The similar behaviors of these three alternating spin 
systems suggest that a ferrimagnet can be considered as a ferromagnet with small 
quantum fluctuations. 



7.2 Low-temperature thermodynamic properties 

We have varied the size of the system from 8 to 20 sites to calculate the thermody- 
namic properties. We impose periodic boundary conditions to minimize finite size 
effects with S^at+i = S^i, so that the number of sites equals the number of bonds. 
We set up the Hamiltonian matrices in the DMRG basis for all allowed S z sectors 
for a ring of 2N sites. We can diagonalize these matrices completely to obtain all 
the eigenvalues in each of the S z sectors. As the number of DMRG basis states 
increases rapidly with increasing m, we retain a smaller number of dominant den- 
sity matrix eigenvectors in the DMRG procedure, i.e., 50 < m < 65, depending on 
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the S z sector as well as the size of the system. We have checked the dependence of 
properties (with m in the range 50 < m < 65) for the system sizes we have studied 
(8 < 2N < 20), and have confirmed that the properties do not vary significantly for 
the temperatures at which they are computed; this is true for all the three systems. 

It may appear surprising that the DMRG technique which essentially targets a 
single state, usually the lowest energy state in a chosen S z sector, should provide 
accurate thermodynamic properties since these properties are governed by the en- 
ergy level spacings and not by the absolute energy of the ground state. However, 
there are two reasons why the DMRG procedure yields reasonable thermodynamic 
properties at low temperatures. Firstly, the projection of the low-lying excited state 
eigenfunctions on the DMRG space which contains the ground state is substantial; 
hence these excited states are well described in the chosen DMRG space. Secondly, 
the low-lying excitations of the full system are often the lowest energy states in dif- 
ferent sectors in the DMRG procedure; hence their energies are quite accurate even 
on an absolute scale. 

The canonical partition function Z for the 2N site ring can be written as 

Z = e*" B(M)l) , (46) 
j 

where the sum is over all the DMRG energy levels of the 2N site system in all the 
S z sectors. Ej and (M)j denote the energy and the ^-component of the total spin 
of the state j, and B is the strength of the magnetic field in units of 1/g/iB (g is the 
gyromagnetic ratio and //# is the Bohr magneton) along the z direction. The field 
induced magnetization < M > is defined as 

<M> = ^ 3 y h . (47) 

The magnetic susceptibility x is related to the fluctuation in magnetization 

X = P[ <M 2 > - <M > 2 ] . (48) 



Similarly, the specific heat Cy is related to the fluctuation in the energy and can be 
written as 

(3 2 [ < E 2 > - < E > 2 ] . (49) 



k 



B 



In the discussion to follow, we present results on the 20-site ring although all 
calculations have been carried out for system sizes from 8 to 20 sites. This is because 
the qualitative behavior of the properties we have studied are similar for all the ring 
sizes in this range for all three systems. 
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The dependence of magnetization on temperature for different magnetic field 
strengths are shown in Fig. 19 for all three systems. At low magnetic fields, the 
magnetization shows a sharp decrease at low temperatures and shows paramagnetic 
behavior at high temperatures. As the field strength is increased, the magnetization 
shows a slower decrease with temperature, and for high field strengths the magne- 
tization shows a broad maximum. This behavior can be understood from the type 
of spin excitations present in these systems. The lowest energy excitation at low 
magnetic fields is to a state with spin s less than Sq- Therefore, the magnetization 
initially decreases at low temperatures. As the field strength is increased, the gap to 
spin states with S > So decreases as the Zeeman coupling to these states is stronger 
than to the states with S < Sq- The critical field strengths at which the magnetiza- 
tion increases with temperature varies from system to system since this corresponds 
to the lowest spin gap of the corresponding system. The behavior of the system at 
even stronger fields turns out to be remarkable. The magnetization in the ground 
state (T = 0) shows an abrupt increase signalling that the ground state at this field 
strength has S z > Sq- The temperature dependence of the magnetization shows a 
broad maximum indicating the presence of states with even higher spin values lying 
above the ground state in the presence of this strong field. In all three cases, the 
ground state at very high field strengths should be ferromagnetic. For the systems 
at such high fields, the magnetization decreases slowly with increase of temperature 
as no other higher spin states lie above the ground state. While we have not studied 
such high field behaviors, we find that the field strength corresponding to switching 
the spin of the ground state sq to s G + 1 is higher for (3/2, 1/2) system compared 
to (3/2, 1) and (1,1/2) systems. The switching field appears to depend on the value 
of \s\ — S2I ■ We see in Fig. 19 that in the (3/2, 1) and (1, 1/2) cases, the ground 
state has switched to the higher spin state at the highest magnetic field strength we 
have studied but in the (3/2, 1/2) case, the ground state has not switched even at 
the field strength indicating that the excitation gap for this system is larger than 
the other two. For the (3/2, 1/2) case, the same situation should occur at very high 
magnetic fields. Thus, we predict that the highest S z is attained in the ground state 
at high magnetic field and this field strength increases with increase in the site spin 
difference \si — s 2 \- 

The dependence of xl '/2N on temperature for different field strengths are shown 
in Fig. 20 for all three systems. For zero field, the zero temperature value of xT is in- 
finite in the thermodynamic limit; for finite rings it is finite and equal to the average 
of the square of the magnetization in the ground state. For the ferrimagnetic ground 
state xT/2N, as T — ■> 0, is given by Sq(Sg + l)/6iV. As the temperature increases, 
this product decreases and shows a minimum before increasing again. For the three 
systems studied here, the minimum occurs at different temperatures depending on 
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the system. For the (3/2, 1) alternating spin system, it is at k B T = (0.8 ± 0.1) J, 
while for the (3/2,1/2) and (1,1/2) cases, it occurs at k B T = (1.0 ± 0.1) J and 
ksT = (0.5 ± 0.1) J respectively. The minimum occurs due to the states with 
S z < Sq getting populated at low temperatures. In the infinite chain limit, these 
states turn out to be the gapless excitations of the system. The subsequent increase 
in the product yT is due to the higher energy-higher spin states being accessed 
with further increase in temperature. This increase is slow in (3/2, 1/2) case, as in 
this system very high spin states are not accessible within the chosen temperature 
range. Experimentally, it has been found in the bimetallic chain compounds that 
the temperature at which the minimum occurs in the yT product depends upon the 
magnitude of the spins s\ and S2 [j39|j . The Ni 11 — Cu 11 bimetallic chain shows a min- 



imum in xT/2N at a temperature corresponding to 55 cm 1 (80K); an independent 
estimate of the exchange constant in this system is 100 cm -1 |JT]|. This is in very 



good agreement with the minimum theoretically found at temperature (0.5 ± 0.1) J 
for the (1,1/2) case. The minimum in yT '/2N vanishes at B — O.lJ/gfis which 
corresponds to about 10T for all three systems. It would be interesting to study 
the magnetic susceptibility of these systems experimentally under such high fields. 
The low-temperature zero-field behavior of our systems can be compared with the 
one-dimensional ferro magnet. In the latter, the spin wave analysis shows that the 
yT product increases as 1/T at low temperatures p2[ . 



In finite but weak fields, the behavior of yT is different. The magnetic field 
opens up a gap and yT falls exponentially to zero for temperatures less than the 
gap in the applied field for all three systems. Even in this case a minimum is found 
at the same temperature as in the zero-field case for the corresponding system, for 
the same reason as discussed in the zero field case. 

In stronger magnetic fields, the behavior of yT from zero temperature up to 
ksT = J m i n (J m in is the temperature at which the minimum in yT is observed) 
is qualitatively different. The minimum in this case vanishes for all three systems. 
In these field strengths, the states with higher S z values are accessed even below 
ksT = J m in- The dependence of yT above k B T = J min at all field strengths is 
the same in all three systems. In even stronger magnetic fields, the initial sharp 
increase is suppressed. At very low temperature, the product yT is nearly zero and 
increases almost linearly with T over the temperature range we have studied. This 
can be attributed to a switch in the ground state at this field strength. The very 
high temperature behavior of yT should be independent of field strength and should 
saturate to the Curie law value corresponding to the mean of magnetic moments 
due to spin-si and spin-S2- 

The temperature dependence of specific heat also shows a marked dependence on 
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the magnetic field at strong fields. This dependence is shown in Fig. 21 for various 
field strengths for all the three systems. In zero and weak magnetic fields, the specific 
heat shows a broad maximum at different temperatures which are specific to the 
system. Interestingly, the temperature at which the specific heat shows a maximum 
closely corresponds to the temperature where the low-field \T has a minimum for 
the corresponding system. For a strong magnetic field (B ~ J), there is a dramatic 
increase in the peak height at about the same temperature corresponding to the 
specific system, although the qualitative dependence is still the same as at low 
magnetic fields in all three cases. This phenomena indicates that the higher energy 
high-spin states are brought to within ksT of the ground state at this magnetic field 
strength for all three cases. 

Studies of thermodynamic properties of the dimerized alternating spin chains 
in these three cases show qualitatively similar trends to that of the corresponding 
uniform systems; this follows from the fact that the low-energy spectrum does not 
change qualitatively upon dimerization. 



8 Magnetization Properties of a Spin Ladder 

As mentioned earlier, a quantum spin system can sometimes exhibit magnetization 
plateaus as a function of an applied magnetic field [37-45] . In this section, we will use 
the finite system DMRG method to study the magnetic properties of a three-legged 
spin-1/2 ladder |]93" |. We consider the Hamiltonian 

3 3 

h = j' s Qin • s a+ i jn + j ^2 ^2 s a , n ■ s ain+ i — u ^2 ^2 s* n , (50) 

a n a =l n a=l n 

where a denotes the chain index, n denotes the rung index, h denotes the magnetic 
field (we have absorbed the gyromagnetic ratio g and the Bohr magneton fiB in the 
definition of h), and J, J' > 0. It is convenient to scale out the parameter J, and 
quote all results in terms of the two dimensionless quantities J'/J and h/J. If the 
length of each chain is L, the total number of sites is N = 3L. Since the total S z 
is a good quantum number, it is more convenient to do the numerical computations 
without including the magnetic-field term in fl5"0|), and then to add the effect of 
the field at the end of the computation. For the ground state properties, we have 
only considered an open boundary condition (OBC) in the rung direction, namely, 
the summation over a in the first term of ([50D runs over 1,2. However, for low- 
temperature properties, we have studied both OBC, as well as a periodic boundary 
condition (PBC) in the rung direction in which we sum over a = 1, 2, 3 in the first 
term. 
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We have done DMRG calculations (using the finite system algorithm [611) with 
open boundary conditions in the chain direction. We have gone up to 120 sites, 
i.e., a chain length of 40. The number of dominant density matrix eigenstates, 
corresponding to the m largest eigenvalues of the density matrix, that we retained 
at each DMRG iteration was m = 80. In fact, we varied the value of m from 60 to 
100, and found that m = 80 gives satisfactory results in terms of agreement with 
exact diagonalization for small systems and good numerical convergence for large 
systems. For inputting the values of the couplings into the numerical programmes, 
it is more convenient to think of the system as a single chain (rather than as three 
chains) with the Hamiltonian 

2 27T? -.a ~ 

H = o J ' E [ 1 ^ cos ( — ) ] S, • S i+1 + J £ S, • S i+3 . (51) 
6 i 6 i 

The system is grown by adding two new sites at each iteration. Note that our method 
of construction ensures that we obtain the three-chain ladder structure after every 
third iteration when the total number of sites becomes a multiple of 6. At various 
system sizes, starting from 48 sites and going up to 120 sites in multiples of 6 sites, 
we computed the energies after doing three finite system iterations; we found that 
the energy converges very well after three iterations. The energy data is used in 
Figs. 22 and 23 below. After reaching 120 sites, we computed the spin correlations 
after doing three finite system iterations. This data is used in Figs. 24 and 25. All 
our numerical results quoted below are for J/ J' = 1/3. We chose this particular 
value of the ratio because there is a particularly broad magnetization plateau at 
m s = 1/2 which can be easily found numerically. 

We now describe the various ground-state properties we have found with OBC 
along the rungs. We looked for magnetization plateaus at m s = 0, 1/2 and 1. For a 
system with N sites, a given value of magnetization per rung, m s , corresponds to a 
sector with total S z equal to M — m s N/3. Using the infinite system algorithm, we 
found the lowest energies E (S Z , N) in the three sectors S z = M + 1, M and M — 1. 
Then we examined the three plots of Eq/NJ vs 1/N and extrapolated the results 
up to the thermodynamic limit iV — > oo. We fitted these plots with the formula 
Eq/NJ = e, + cti/N + bi/N 2 , where the label i = 1,2,3 denotes the S z sectors 
M+l, M and M — 1. In the thermodynamic limit, the values of the three intercepts 
Ci should match since those are just the energy per site for the three states whose 
S z, s differ by only 1. However, the three slopes are not equal in general. We now 
show that there is a magnetization plateau if a\ + 03 — 2a2 has a nonzero value. Since 
the three energies Eq are computed without including the magnetic field term, the 
upper critical field h c+ where the states with S z = M + 1 and M become degenerate 
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is given by 

h c+ {N) = E {M+1,N) - E {M,N) . (52) 

Similarly, the lower critical field /i c _ where the states with S z = M and M — 1 
become degenerate is given by 

K4N) = E (M,N) - E (M-1,N). (53) 

We therefore have a finite interval Ah(N) = h c+ (N) — h c ^(N) in which the lowest 
energy state with S z = M is the ground state of the system with N sites in the 
presence of a field h. If this interval has a nonzero limit as N — > oo, we have a 
magnetization plateau. Thus, in the thermodynamic limit, the plateau width Ah/ J 
is equal to ai + a 3 — 2a 2 . 

We will now quote our numerical results for J/ J' = 1/3. For a rung magneti- 
zation of m s = 1/2, i.e., M = N/6, we found the three slopes to be equal to 
3.77, —0.02 and —1.93; see Fig. 22. This gives the upper and lower critical fields to 
be 

h c+ /J = ai — a 2 = 3.79 , h c _/J = a 2 — a 3 = 1.91 , 

Ah/ J = (h c+ - h c J)/J = 1.88 . (54) 

This is a sizeable plateau width. For a rung magnetization of m s = 1, we found the 
^ to be equal to 4.97, —0.24 and —5.43. Thus the upper and lower critical fields are 

h c+ /J = 5.21 , h c _/J = 5.19 , Ah/ J = 0.02 . (55) 

Finally, for a rung magnetization of m s = 0, we need the energies of states with 
M = and M — ±1. Since the last two states must have the same energy, we have 
ai = a 3 and it is sufficient to plot only E (0, N) and E (l, N) vs \/N . We found a± 
and a 2 to be equal to 0.39 and 0.34. This gives the upper and lower fields to be 

h c+ /J = 0.05, h c ^/J = -0.05, Ah/ J = 0.10. (56) 

The plateau widths given in (|55D and (|56|) are rather small. In Fig. 23, we indicate 
the plateau widths Ah(N)/J as a function of 1/N for m s = 1/2, and 1. 

Next, we computed various two-spin correlations for the 120-site system. These 
are denoted by (S^S^) and (S^S^). For the zz correlations, it is convenient to 
subtract the product of the two separate spin densities. At m s = 1/2, we found that 
all these correlations decay very rapidly to zero as the rung separation \l — n\ grows. 
In fact, the fall offs were so fast that we were unable to compute sensible correlation 
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lengths. The correlation lengths are of the order of one or two rungs as can be seen 
in Fig. 24. 

On the other hand, for the state at m s = 0, we found that all the correlations 
decay quite slowly. The decays are consistent with power law fall offs of the form 
A(— l)''~ n '/|Z — n\ v . It is difficult to find rj very accurately since the maximum value 
of | Z — nj is only 20; this is because we fixed one site to be in the middle of the chain 
(to minimize edge effects), and the maximum chain length is 40 for our DMRG 
calculations. For m s = 0, the exponent r\ for all the correlations was found to be 
around 1. There was no difference in the behaviors of the zz and H — correlations 
since this was an isotropic system; m s = is the ground state if the magnetic field 
is zero. 

For the state at m s = 1 (which is the ground state only for a substantial value 
of the magnetic field), we found that the H — correlations again decay quite slowly 
consistent with a power law. The exponents r] for the different H — correlations 
varied from 0.61 to 0.70 with an average value of 0.66; see Fig. 25 for an example. 

We now describe some low-temperature thermodynamic properties of the three- 
chain system obtained using DMRG. Although DMRG is normally expected to be 
most accurate for targeting the lowest states in different S z sectors, earlier stud- 
ies of mixed spin chains have shown that DMRG is quite reliable for computing 
low-temperature properties also |78| . There are two reasons for this; the low-lying 



excited states generally have a large projection onto the space of DMRG states which 
contains the ground state, and the low-lying excitations in one sector are usually 
the lowest states in nearby S z sectors. 

We first checked that for systems with 12 sites, the results obtained using DMRG 
agree well with those obtained by exact diagonalization. We then used DMRG to 
study the magnetization, susceptibility and specific heat of 36-site systems with both 
OBC and PBC along the rungs. We computed these quantities using the techniques 
described in Eqs. ( ^6] - |49|) . The plots of magnetization vs magnetic field for various 
temperatures are shown in Fig. 26 for OBC along the rungs. The temperature T 
is measured in units of J/k B . We see that the plateau at m s = 1/2 disappears 
quite rapidly as we increase the temperature. The plateau has almost disappeared 
at T = 0.4 which is substantially lower than the width Ah /J = 1.88. The mag- 
netic susceptibility is (exponentially) small at low temperatures in the region of the 
plateau because the magnetic excitations there are separated from the ground state 
by a gap. Similar results for the magnetization and susceptibility are found for the 
case of PBC along the rungs. 

However, the specific heats at the m s = 1/2 plateau demonstrate an interesting 
difference between OBC and PBC along the rungs. While it is very small at low 
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temperatures for OBC (see Fig. 27), it is not small for PBC, although it shows a 
plateau in the same range of magnetic fields as the magnetization itself. This obser- 
vation strongly suggests that the system with PBC along the rungs has nonmagnetic 
excitations which do not contribute to the magnetization or susceptibility, but do 
contribute to the specific heat. Fig. 28 gives a more direct comparison between 
OBC and PBC along the rungs. Although these nonmagnetic excitations were stud- 
ied by previous authors |58], [94]] , we believe that our specific heat plots prove 
their existence most physically. To show these excitations even more explicitly, we 
present in Fig. 29 all the energy levels for a 12-site chain in the sector S z = 2 
(i.e., m s = 1/2) using exact diagonalization. It is clear that the ground state is 
well separated from the excited states for OBC, but it is at the bottom of a band 
of excitations for PBC; these excitations are nonmagnetic since they have the same 
value of S z as the ground state. 

We summarize our results for a three-chain spin- 1/2 ladder with a large ratio 
of interchain coupling to intrachain coupling. There is a wide plateau with rung 
magnetization given by m s = 1/2 for both OBC and PBC along the rungs. For the 
case of OBC, the two-spin correlations are extremely short-ranged, and the magnetic 
susceptibility and specific heat are very small at low temperature in the plateau. All 
these are consistent with the large magnetic gap. At other values of m, the two-spin 
correlations fall off as power laws. For the case of PBC, the magnetic susceptibility 
is again very small at low temperature in the plateau. However the specific heat goes 
to zero much more slowly which dramatically shows the presence of nonmagnetic 
excitations. 

To summarize this review, we have discussed a variety of numerical and analytical 
methods for studying spin clusters and quasi-one- dimensional spin systems. The 
methods discussed for spin clusters are directly relevant to areas of current interest 
such as quantum tunneling in the presence of time- dependent magnetic fields [17- 
20,95]. One-dimensional spin systems sometimes have unusual properties (such as 
a disordered ground state with an excitation gap above it) which are not observed 
in similar models in higher dimensions. The techniques described here, particularly 
the DMRG method, are well-suited for studying extended systems in one dimension. 
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Table 1: Low-Lying states of Mn 12 Ac, relative to the ground state for the param- 
eter in question. Entries in parenthesis in cases A, B and C correspond to the 



effective Hamiltonian results of Sessoli et al. f22|. Case D corresponds to the 
parameters suggested by Chudnovsky [|T7|]. The parameters corresponding to dif- 
ferent cases are: case (A) J 1= 225K, J 2 =90K, J 3 =90K, J 4 =0K; case (B) J 1= 225K, 
J 2 =90K, J 3 =93.8K, J 4 =0K; case (C) J 1= 225K, J 2 =90K, J 3 =86.2K, J 4 =0K; case 
(D) J 1= 215K, J 2 =85K, J 3 =-85K, J 4 =-45K; case (E) J 1= 215K, J 2 =85K, J 3 =85K, 
J 4 =-64.5K. All the energies are in K. 
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Table 2 : Energies (in unit of K) of few low lying states in Fe§. The exchange 
constants corresponding to the various cases are: case (1) Ji = 150K, J 2 = 25K, J 3 
= 30K , J 4 = 50K; case (2) Ji = 180K, J 2 = 153K, J 3 = 22.5K , J 4 = 52.5K; case 
(3) Ji = 195K, J 2 = 30K, J 3 = 52.5K , J 4 = 22.5K. All the energies are in K. 
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Figure Captions 



1. A schematic representation of the VB diagrams for eight spin-1/2 objects. 
States |1 > and |2 > are legal VB diagrams and |3 > is an illegal VB diagram. 

2. VB diagram for spin-1/2 objects with total spin 1/2 (|1 >) and total spin 1 
(|2 >). Their bit representations (0 and 1) as well as the unique integer 
representing them are shown. P and P' are the phantom sites. |3 > is a singlet 
VB diagram corresponding to two spin-1, a spin-5/2 and a spin-3/2 object. 

3. Effect of operation by the operator (S« ■ Sj — 1/4) on a state with a singlet line 
between sites i and j and on a state with sites i and j singlet paired with two 
different sites i' and f. 

4. A schematic diagram of the exchange interactions between the Mn ions in the 
Mn 12 Ac molecule. 

5. Spin densities in the ground state (S=10, M s =10) of Mn 12 Ac for parameter 
values Ji=215K, J 2 =85K, J 3 =85K and J 4 =-64.5K. 

6. A schematic diagram of the exchange interactions between the Fe ions in the 
Fe 8 molecule. 

7. Spin density in the ground state (S=10,M S =10) of Fe 8 for three different pa- 
rameter values: (a) J 4 = 150K, J 2 = 25K, J 3 = 30K, J 4 = 50K, (b) J x = 180K, 
J 2 = 153K, J 3 = 22.5K, J 4 = 52.5K, and (c) J 4 = 195K, J 2 = 30K, J 3 = 52.5K, 
J 4 = 22.5K. 

8. Schematic picture of the frustrated and dimerized spin chain. 

9. Classical phase diagram of the spin chain in the J 2 — 5 plane. 

10. 'Phase' diagram for the spin-1 chain in the J 2 — 5 plane. 

11. Doubly degenerate ground states (a) and (b) of the J — J 2 chain (see Fig. 8 for 
5 = 0) at J 2 = J/2. The solid line between sites i and j represents a singlet, 

[I u s > -I uh>]/V2. 

12. 'Phase' diagram for the spin-1/2 chain in the J 2 — 5 plane. 
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13. Gap A vs J for coupled spin chains (5 = 1). Spin- 1/2 and spin-1 data are 
indicated by crosses and circles respectively. 

14. Expectation values of the z-components of the two spins vs the unit cell index 
n for an alternating spin chain. The upper and lower points are for the spin-si 
and the spin-s 2 sites respectively. 

15. Subtracted two-spin correlation functions as a function of distance between the 
two spins, (a) spin-Si spin-Sx correlations, (b) spin-s 2 spin-s 2 correlations, and 
(c) spin-Sx spin-s 2 correlations. In each figure, squares correspond to (3/2, 1), 
circles to (3/2, 1/2) and triangles to (1, 1/2) systems. 

16. Energy difference (units of J) between the ground state and the lowest energy 
state with spin S = So — 1 as a function of inverse system size. Sq is the total 
spin of the ground state. 

17. Excitation gap (units of J) from the ground state (spin S = So) to the state 
with spin S = Sg + 1 as a function of the inverse system size. 

18. Excitation gap (units of J) to the state with spin S — Sq + 1 from the ground 
state (S = Sq) as a function of 5 for the dimerized alternating chain. The 
exponent is 1.0 ± 0.01 for all three systems. 

19. Plot of magnetization per site as a function of temperature T for four different 
values of the magnetic fields B. Squares are for B = 0.1J/g/j, B , circles for 
B = 0.5J/g/i B , triangles for B = J/g/iB and diamonds for B = 2J/gfi B . 

20. yT ( defined in the text) per site as a function of temperature T for various 
magnetic fields B. Zero field results are shown by squares, B = 0.01J/gfj, B by 
circles, B = O.U/g/iB by triangles and B = J/g/iB by diamonds. 

21. Specific heat per site as a function of temperature T for four different values 
of magnetic fields B. Zero field data are shown by squares, B = 0.01J/gfi B 
by circles, B = 0.1J/g/i B by triangles and B = J j g\i B by diamonds. 

22. The energy/site in units of J vs l/N at the m s = 1/2 plateau, for the three- 
chain ladder with J/ J' = 1/3. The curves indicate quadratic fits for (a) 
E (M + 1, N), (b) E (M, N), and (c) E (M — 1, TV). 

23. Plateau widths vs l/N for (a) m s = 1/2, (b) m s = 0, and (c) m s = 1. 

24. Correlation function < S^S^n > at the m s = 1/2 plateau for J/ J' = 1/3. 
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25. Correlation functions < Si l S 2n > in the m s — 1 state for J/ J' = 1/3. 

26. Magnetization vs magnetic field at six different temperatures, for a 36-site 
system with OBC along the rungs and J/ J' = 1/3. 

27. Specific heat in units of ks vs magnetic field at six different temperatures, for 
a 36-site system with OBC along the rungs and J/ J' = 1/3. 

28. Comparisons of specific heat and susceptibility vs the magnetic field for a 
36-site systems with OBC and PBC along the rungs. 

29. Comparison of the energy spectra in units of J of the 12-site system with OBC 
and PBC along the rungs. The energies in the S z = 2 sector are shown for 
J/ J' = 1/3. 
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